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, Abstract 

In a recent series of papers, the class of energy-conserving Runge-Kutta methods 
. named Hamiltonian BVMs (HBVMs) has been defined and studied. Such methods 

\ have been further generahzed for the efficient solution of general conservative prob- 

r~| ' lems. In this paper we derive a further extension, which we name Enhanced HBVMs 

"j^ ■ (EHBVMs), more tailored for Hamiltonian problems, allowing for the conservation of 

multiple invariants of the continuous dynamical system. The analysis of the methods 
is fully carried out and some numerical tests are reported, in order to confirm the 
^ ^ theoretical achievements. 
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^ : 1 Introduction 

Hamiltonian problems arise in many field of application, ranging from the nano-scale of 
^ ! molecular dynamics to the macro-scale of celestial mechanics. Such problems are in the 
following form: 

y' = JVHiy), y{0) = y^ e M^-, (1) 
where the state vector is often partitioned as 



y = \ p ] ^ (i^p^ 



nm 



Moreover, 



J={ ]=-J^ = -J-\ (2) 
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and H{y) = H{q,p) is the Hamiltonian function defining the problem. From ([T]) and it 
is straightforward to derive that H{y{t)) = H{yo) for t > 0, since 

^^Hivit)) = VH{y{t)fy'{t) = V H{y{t)f JV H{y{t)) = 0, 

due to the fact that J is skew- symmetric. For isolated mechanical systems, the Hamiltonian 
has the physical meaning of the total energy of the system and, therefore, it is of interest 
to derive methods which are able to preserve this property in the discrete solution. For the 
continuous problem, it can be seen that the symplecticity of the map implies the property of 
energy conservation of the given system, so that a relevant line of investigation, concerning 
the efficient numerical solution of such problems, has been that of devising symplectic meth- 
ods, namely methods for which the discrete map inherits the property of symplecticity (see, 
e.g., [291 EQ]). In particular, in [29] the existence of infinitely many symplectic Runge- 
Kutta methods was proved, and an algebraic criterion for symplectic Runge-Kutta methods 
is provided. 

Nevertheless, unless the continuous case, in the discrete setting the symplecticity of the 
map does not imply energy-conservation (see also [H]), so that a different line of investigation 
has been that of looking for energy-conserving methods. One of the first approaches along 
this line is represented by discrete gradient methods fl8\ |27] . which are based upon the 
definition of a discrete counterpart of the gradient operator, so that energy conservation 
for the numerical solution is guaranteed at each step and for any choice of the integration 
step-size. A different approach is based on the concept of time finite element methods, which 
has led to the definition of energy-conserving Runge-Kutta methods lU [2l |3T1 |32], based 
on a local Galerkin approximations for the equation. A partially related approach is given 
by discrete line integral methods [211 [23 [26j, where the key idea is to exploit the relation 
between the method itself and the discrete line integral, i.e., the discrete counterpart of 
the line integral in conservative vector fields. This, in turn, allows exact conservation for 
polynomial Hamiltonians of arbitrarily high-degree, resulting in the class of methods later 
named Hamiltonian Boundary Value Methods (HBVMs), which have been developed in a 
series of papers [Bl [HI [IHl [Zl [HI [121 [131 [H [S] (we refer to [B] for a systematic presentation of 
this approach). Another approach, strictly related to the latter one, is given by the averaged 
vector field method [2H] and its generalizations [TH] , which have been also analysed in the 
framework of B-series [161 [22] (i.e., methods admitting a Taylor expansion with respect to 
the step-size). 

For sake of completeness, we also mention that attempts aimed to obtain methods that, 
in a weaker sense, have both the property of symplecticity and energy-conservation have been 
also considered (see, e.g., [2S1 [IH [33]). 

Sometimes, the dynamical system defined by ([1]) has additional invariants, besides the 
Hamiltonian. It is therefore interesting to devise methods which are able to preserve all 
of them in the discrete solution. The approach based on the discrete line integrals, which 
HBVMs rely on, has been then used to cope with this problem, leading to the class of line 
integral methods which are able to preserve any number of invariants for general conservative 
problems [3] (see also [6]). In this paper, we consider a different generalization of HBVMs, 
still based on the concept of discrete line integral, which is able to provide multiple invariant- 
conserving methods when the problem is in the form ([1]). With this premise, the paper is 
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organized as follows: in Section [2] we recall the basic facts about HBVMs; in Section E] we 
define their multiple invariants- conserving extension; in Section |4] we provide numerical tests 
for the presented methods; finally, in Section [5] we give some conclusions. 



2 HBVMs 



Let us consider a polynomial approximation to the solution of ([T]) over the interval [0,h], in 
the form 



s-l 



(3) 



a'{ch) = Y,P,{c)j,ia), cG [0,1], 
i=o 

where {Pj}j>o is the family of Legendre polynomials, shifted and scaled in order to be 
orthonormal on the interval [0, 1], 

»i 



deg Pj = j, 



Pj{x)Pj{x)dx = 6, 



Vz,j>0. 



(4) 



By imposing the initial condition cr(0) = i/q, and setting yi = a{h) ^ y{h), the coefficients 
7j(cr) are determined by imposing the conservation of energy sX t = h. This, in turn, is 
gained by imposing that the following integral vanishes: 



= H{yi) - H{y^) = H{a{h)) - H{a{0)) = / VH{a{t)fa'{t)dt 

Jo 

= h [ VH{a{Th)fa\Th)dT. 
Jo 



(5) 



By taking into account of Q, one then requires [TO] : 



s-l 

E 

j=0 

which holds true, provided that 



Pj{T)VH{a{Th))di 



7,(a)=0, 



(6) 



lA^)=V,J f P,{T)VH{a{Th))dT, 
Jo 



J = 0,...,s-1, (7) 

where tjq, . . . , r]s_i are arbitrary constants. HBVMs are then obtained by setting 

= 1) J = 0, . . . , s - 1, 
resulting in an approximation of order 2s to y{h) [TUl fT2] : 

a{h)-y{h) = 0{h'^^'). 
In particular, by considering the orthonormality of the polynomial basis, one obtains that 



yi = a{h) = yo + h-fo{a) = yo + j JVH{a{t))dt. 
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This latter expression, clearly show that this polynomial approximation generalizes that 
defined in [28]. However, the resulting polynomial approximation, given by 

aich)=yo + hJ2 Pii^)'^^lji^)^ cG [0,1], (8) 
i=o 

provides an effective numerical method, only when the integrals appearing in ([7]) are conve- 
niently approximated by means of a quadrature formula. If this latter formula is defined at 
the k Gauss- Legendre points in [0, 1], 



(i.e., Pk(ci) = 0, i = 1, 



< ci < . . . < Cfc < 1, 
k) and corresponding quadrature weights 
6i, . . . > 0, 



(9) 



(10) 



one then obtains a HBVM(fc, s) method which can be cast as a fc-step Runge-Kutta method, 
with abscissae (Q, weights ffTOj) . and Butcher matrix given by 



:iii 



where 



V, = (P,-i(q)), Is 



Pj_i(x)dx 1 G 



nkxs 



Q = diag(6i. 



The corresponding polynomial approximation is then given by [TOj [T2 



u{ch) 



S-l „c 

+ h'^ Pj{x)dx 



j=0 
s-l 



J2biPjici)JVH{ 



= Vo + h^ I Pj{x)dx%, cG [0,1], 

j=0 



where 



:i2) 



(13) 



Ui = u{cih), £=l,...,k, 

are nothing but the stages of the Runge-Kutta method. It can be proved that [H [lOl [12], for 
all k>s, a HBVM(fc, s) method: 

• has order 2s; 

• is symmetric; 

• when A; = s it reduces to the s-stage Gauss- Legendre method; 

• is energy-preserving for all polynomial Hamiltonians of degree not larger than 2k/ s. 
Differently, the error in the Hamiltonian is 0{h'^''~^^), provided that H is suitably reg- 
ular. 
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From the last point, a practical conservation of the Hamiltonian follows, also considering 
that the computational complexity of the method is s, independently of k. Indeed, by 
reformulating the discrete problem generated by the method in terms of the s unknown 
coefficients {7^} appearing in ( !T2|) . one obtains the discrete problem [IT] 



7 = Vjn ® / JVH {e®yo + hi, ® I ^) 



where 



/ 1 \ 



7 



/ 7o ^ 
V 7. 1 / 



which has (block)-size s, independently of k. 

It is worth mentioning that, because of the existing relations between the integrals of the 
polynomials {Pj} and the polynomials themselves, matrix f lTT]) can be also written as 



(14) 



where 



2 
^1 







V 



... 



with 



By considering that, for k = s: 

Vs+i ={Vs ) , 



vjn = V. 



1 

S 1 



one then sees that (fT4|) can be regarded as a generalization of the ly-transformation for 
collocation methods, as defined by Hairer and Wanner [211 page 79] . 



3 Mutiple invariants preserving HBVMs 

We now use again the approach based on line integrals, to define a multiple-invariant con- 
serving version of HBVM(A;, s) methods. Though the basic idea is similar to that used in 
|3], nevertheless, the obtained methods are definitely different from those described in that 
reference. 

In more detail, we now use the fact that energy conservation is gained, in ([7]), whichever 
r]j is. Assume then that 

L : M^™ ^ M'^ (15) 

is a set of (functionally independent) smooth invariants for the dynamical system ([1]). Con- 
sequently, one has 

VL{yfJVH{y) = G M^ Vy, (16) 
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where \/L(y)'^ is the Jacobian matrix of L. We will now extend the approach described in the 
previous section, for imposing their conservation, assuming that u < s. For sake of simplicity, 
we shall at first define a polynomial approximation a G 11^ at a continuous level (i.e., similar 
to ([S])-®), then passing to define a fully discrete approximation u Ellg. Clearly, by setting 
a in the form ([3]), we gain energy-conservation by repeating similar steps as done up to ([HD. 
The difference, in such obtained by setting 



Vj 



Vj 



1, j = 0, . . . ,s - - 1, 



[1 - 



(17) 



j = s-z/, 



in order to obtain the conservation of the u additional invariants ( !T5|) - (fT6|) . even though, in 
principle, any subset of u coefficients {rjj} could be used for this purpose. By setting, as 
before, 

yi = (^(h) ^ y{h) 

the new approximation, from ([3]), ([7]), and f|T7|) one then obtains that: 

= L{y,) - L{yo) = L{a{h)) - L{a{0)) = / VL{a{t)fa'{t)dt 

Jo 



»i s—l p „i 

h / VL((T(r/i))V(r/i)dr = hV] / P,(r) VL((T(r/i))dr 
Jo Uo 



7j(^)dr 



s-l 



_j=0 j=s-v 

where (see ([7]) and (fTTll ). for all j > 0: 



cr 



I Pj{T)VL{a{Th)y\T Gl 
Jo 

[ Pj{T)JVH{a{Th))dT G 
Jo 



2mxu 



X,2m 



Consequently, energy-conservation is "for free" and, moreover, the invariant conservation is 
gained provided that 



s-l 



s-l 



J2 h'(^-'-^^a,<j>,{ar^,{a) 



J = S-U 



j=0 



By defining the matrix 

r(a) = [ /i2(-i)0,_,(a)^7s-.(a), 
and the vectors 



<i_i((T)^7s_i(a) 



] G R"' 



(19) 



(20) 



, b(a) = 5^0,(a)%-(a) GM^ 



j=0 



equation fll9l) can be recast in vector form as 

T{a)a = h{a). (21) 

The following results then hold true. 

Lemma 1 Let ip : [0, h] — > V , with V a vector space, admit a Taylor expansion at 0. Then, 
for all j > 0: 

-I 

Pj{T)^lj{Th)dT = 0{y). 
Proof. By taking into account (jl]), one obtains: 

Pj{T)i){Th)dT 



"^0 n>0 ■ n>0 ■ 



Lemma 2 The right-hand side of problem ([2]j can he expanded as 

JVH{y{ch)) = 5^P,(c)7,(l/), c G [0, 1], 

i>o 

where 7(1/) is defined according to ^E), provided that H is suitably regular. 
Proof. See tlSj.D 

Lemma 3 With reference to f21\) . one has: h{a) = 0{h'^^). 
Proof. From (fT6l) and ( ITSl) one obtains: 

^0,(af7,(a) = O. 

Consequently, by virtue of Lemma [H 

j=0 j>s 

Lemma 4 Matrix r(cr) in l[21\) has 0(/i^*^^) entries. 

Proof. The proof follows immediately from ( l20l) and Lemma [H □ 

In order to simplify the subsequent arguments, we make the following assumption on 
matrix r(cr)j3 



^Actually, it suffices that system (|21l) is consistent, but the arguments become more involved. 



Assumption 1 Matrix T{a) is nonsingular. 

The following result then follows easily from Lemmas [3] and |H 

Theorem 1 Under AssumptionUl the vector ex. in l[21\) has 0{h'^) entries. 

Remark 1 We observe that, in order for TheoremU\ to hold, it is necessary that 

V <s. (22) 

In fact, when v = s, from / fTPj) one obtains that the products h'^^'^~~^~^^aj are all equal to 1 
and then, from p7| j it follows that r^j = 0, j = 0, . . . , s — 1. Consequently, in the sequel we 
shall assume that 

Vo = 1. (23) 

We can now state the following result. 

Corollary 1 Under Assumption Ql the method conserves all the invariants. Moreover, 
a{h)-y{h) = 0{h^'+^). 

Proof. The first part of the proof derives from the very definition of the method. The 
second part of the proof strictly follows the technique used in [13]. Let then y(t;co,z) be 
the solution of problem ([1]) satisfying the initial condition y{oj) = z. Moreover, let $(t, r) 
be the corresponding fundamental matrix solution of the associated variational problem. 
Consequently, from Lemmas [T] and 121 Theorem [H and from (fTTll . one obtains: 



aih)-yih) = yih;h,aih))-yih;0,a{0)) 



-y{h;t,a{t))dt 



^ ^ d d 

y{h; u, a{t))\^^^ + — y{h; t, z) \,^. a\t) 



dtu'"' ' ' ' ''"^=* dz 
$(/i, t) [-JVH{a{t)) + a{t)] dt 



dt 



h / ^{h,Th) 
Jo 

-h [ <^{h,Th) 
Jo 



5^P,(r)7,(a) + 5^P,(r)7,(a) dr 

j>0 j=0 
,j>s j=s~v 



dr 



"E 



]=S-V 



Pj{T)(^{h,Th)dT 



=0{h3) 



Pj{T)^{h,Th)dT 



=0{ho) 



=0{hi) 

7i(^) - 



=0{h^) 

=0(/i^) 



2s+l- 
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3.1 Discretization and EHBVM(A;, s) methods 

As is clear, the polynomial approximation cr G 11^ defined above doesn't yet define a numerical 
method: this will be obtained once the integrals in (fT8|) are approximated by means of a 
suitable quadrature formula. As in the case of HBVM(A;, s) methods, for this purpose we 
choose the abscissae placed at the k Gauss- Legendre points in [0,1], and the corresponding 
weights fllOp . In so doing, we obtain a new polynomial approximation, say m G H^, defined 
by replacing the integrals with the given quadrature formula, having order 2k. By setting ui 
formally as in f|T3l) . for < j < s — 1 one then obtains: 



7j 



k 

^htPj{cf)JVH{u^) = 7,(m) - A,(/i), 



(24) 



in place of f|T8l) where, setting 

and assuming L and H suitably regular 



2k 



0, 



if L G n 



and 



0, 



otherwise, 

if HeU^, 
otherwise. 



(25) 



(26) 



(27) 



Setting by % and aj, respectively, the discrete approximations to ( fTTj) . the new polynomial 
approximation is then given by 



u{ch) = Ho + h'S^ / Pj{x)dx fjj^j 



(2J 



i=o 



yo + h 



rP,(x)dx7,- rPA^)dxh'^^~''^^a,^, 



with the scalars cnj satisfying the equation (compare with (fT9l) ): 



s-1 



j=0 



j=s-v 

Similarly as previously done in fl20l) - fl2Tl) . by defining the matrix 
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(29) 



(30) 



and the vectors 



OL 



equation ( l29l) can be recast in vector form as 



j=0 



Vol = h. 



(31) 



Since the number of the additional invariants ( |T6l) has to satisfy ( 122|) (and, then, ( l23i) holds 
true), similarly as in the case of HBVM(/c, s), the new approximation is given by 



yi = u{h) = yo + h% = yo + h^b(,JV H{ui). 



(32) 



Definition 1 We shall denote by EHBVM{k, s) (Enanced HBVM(k, s) ) the methods defined 
by m-WW- 

The following results then easily follow, providing a discrete counterpart of Theorem [H 

Theorem 2 Under Assumption^ matrix f is nonsingular for all sufficiently small step-sizes 
h, and the vector a. has 0{h^) entries. 

We can now state the following results, concerning the order of accuracy of the discrete 
solution, as well as of the invariants, provided by EHBVM(A;, s) methods. 

Theorem 3 Assume that the Hamiltonian function defining problem (QP is a polynomial of 
degree less than or equal to fi as defined in / f^) . Then, EHBVM{k, s) method is energy- 
conserving. Differently, for all general and suitably regular H, one obtains 



H{y,) - H{yo) = 0{h 

Proof. One has: 

H{y,)-H{yo) = H{u{h)) - H{u{0)) = 



2k+l\ 



VH{u{t)fu{t)dt 



„1 „i s-1 

h / VH{u{Th)fu{Th)dT = h Vif(u(r/i))^ V P,(r)%7jdr 

Jo Jo j^Q 



s-1 
j=0 



VH{u{Th))P,j{T)dT 



s-1 



h^Vjlj{ufJlj = (*) 

j=0 



The first part of the proof easily follows from the fact that, if G 11^, then 

7i = 7j(^i), j = 0,...,s-l. 
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so that (*) = 0, since J is skew-symmetric. In general, assuming that H is suitably regular, 
one has (see (1211) and (E 



Using similar arguments, by means of fl26l) it is possible to prove the following result. 



Theorem 4 Assume that the invariants / fig) of problem are polynomials of degree less 
than or equal to fi as defined in [2B^) . Then, EHBVM{k, s) method is invariants- conserving. 
For all general and suitably regular L, one obtains 

L{y{)~L{y,)=0{h^'^'). 

Next result concerns the order of accuracy of the numerical solution. 

Theorem 5 Assuming that both H and L are suitably regular, for all k > s, the numerical 
solution generated by a EHBVM{k, s) method satisfies 



y,~y{h) = 0{h'^^'). 



That is, the method has order 2s. 



Proof. The proof proceeds in a similar way as that of Corollary (U and by using the same 
notation. One has then: 



yi-y{h) = y{h]h,u{h)) - y{h,0,u{0)) 



dt 



y{h-t,u{t))dt 



dt 



— y{h; u, + ^ yih; t, z) \^^^^^^ u'{t) 

<^{h, t) [-JVH{u{t)) + u{t)] dt 



s-l 

, i>o j=o 

s-l 

s-l s-l 
j=0 j>s j=s-y 



h / $(/i,r/i) 
Jo 

h ! <^{h,Th) 
Jo 

h ! ^{h,Th) 
Jo 



dr 



dr 



s-l 
j=0 



P,(r)$(/i,r/i)dr 



=0(/ii) 



j>s 



Pj{T)<i>{h,rh)dT 



=Oih3) 



=0{h3) 

7jH - 
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s-1 
j=s-u 



Pj{T)^{h,Th)dT 



0{h 



O{ho) 
2s+l\ I ^/T,2s+1 



=0{h^) 

=0{h3) 



For sake of completeness, we also mention the following result, whose proof is straight- 
forward (see, e.g., [HIE]). 

Theorem 6 Provided that the abscissae ^ are symmetrically distributed in the interval 
[0,1], an EHBVM{k, s) method is symmetric. 

3.2 Runge-Kutta type formulation of EHBVM(/c, s) methods 

Though the method fl28|) - fl32|) is not strictly a Runge-Kutta method, nevertheless, it admits 
a Runge-Kutta type formulation which is quite useful to represent it. In more details, we 
already saw that a HBVM(/c, s) methods is a fc-stage Runge-Kutta method defined by the 
following Butcher tableau (see (|TT1) ) 



where, as usual, c is the vector of the abscissae and b is the vector of the weights. Moreover, 
we recall that the only formal difference between a HBVM(/c, s) method and an EHBVM(/c, s) 
method stems from the coefficients rji, ■ ■ ■ ,'fis-i which may assume values different from 1 
(indeed, 170 = 1, as stated in fl23|) ). Consequently, by introducing the diagonal matrix 

= diag(l,r)i, . . . ,r),_i), 

one obtains the following Runge-Kutta type formulation of an EHBVM(/c, s) method: 



As an example, HBVM(2,2) is the usual 2-stage Gauss method, whereas EHBVM(2,2) is 
given by 



1 v/3 

2 6 

1 _L \/3 

2 6 



V3 
12 



As expected, when 171 = 1 one retrieves the usual 2-stage Gauss method. 



4 Numerical tests 

We here show a few numerical tests, based on the well known Kepler problem 
by the Hamiltonian 



H{q,p) 



1 



q,p e 



13], defined 
(33) 
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which has a periodic orbit of period 27r, when the initial condition is chosen as 

(g(f,Po)= ( 0, 0, y^), ee[0,l), 

which is given by an ellypse of eccentricity e, in the g-plane. This problem admits two addi- 
tional invariants of motion, besides the Hamiltonian ( l33l) . given by the angular momentum 

L^{q,p) = q^J,p, ■J2=(^\ ly (34) 

and the Laplace- Rung e-Lenz (LRL) vector, resulting in the following conserved quantity: 

Uq,p) = ie^p)L,iq,p)--fl, (35) 

where, as usual, 61,62 G are the two unit vectors. 

We solve this problem, considering an eccentricity e = 0.6, by using the following methods: 

• the symplectic 3-stage Gauss method; 

• the (practically) energy-conserving HBVM(12,3) method; 

• the EHBVM(12,3) method where it is imposed only the (practical) conservation of the 
angular momentum (IMIl . besides the Hamiltonian (!33|) : 

• the EHBVM(12,3) method where it is imposed both the (practical) conservation of the 
angular momentum (IMl) and of the LRL vector (l35!l . besides the Hamiltonian (!33l) . 

In Table [1] we list the measured errors after 10 periods, thus confirming that all methods 
are sixth-order. Moreover, in Table [2] we list the maximum norm for the vector a defined in 
(I3T1) . over the same interval, for the EHBVM(12,3) method: 

• by imposing only the conservation of the angular momentum, besides the Hamiltonian. 
Here, a^^^ = max„^-^ z ||a„||oo; 

• by imposing both the conservation of the angular momentum and of the LRL vector, 

(2,'] 

besides the Hamiltonian. Here, a, = max„_i t an Lo- 

The obtained results confirm that the entries of the vector a. are actually 0{h'^), as predicted. 
For sake of completeness, in Figured] we plot the two components of the vector a. in the second 
case, when a step-size /i = 7r/30 is used: their periodic behaviour, in accordance with that of 
the solution, is clearly evident. 

At last, concerning the conservation of the invariants, by using a constant step-size h = 
0.1, we have solved the problem over the interval [0, 10^], obtaining the following results: 

• concerning the conservation of the Hamiltonian fl33|) . all methods are (practically) 
energy-conserving, except the symplectic 3-stage Gauss method. However, the Hamil- 
tonian error turns out to be bounded, as expected, as confirmed by the plot in Figure [21 
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Table 1: Error in the numerical solution for the 3-stage Gauss method (Eg), HBVM(12,3) 
{Eh), and EHBVM(12,3) with only angular momentum (-E^"*), and with both angular mo- 

(2) 

mentum and LRL vector {E^ ) conserved. 



h 


Eg 


ratio 


Eh 


ratio 


Ee 


ratio 


^E 


ratio 


7r/30 


1.942e-03 




4.587e-05 




1.017e-05 




1.928e-05 




7r/60 


2.817e-05 


68.9 


7.375e-07 


62.2 


1.644e-07 


61.9 


3.052e-07 


63.2 


7r/120 


4.346e-07 


64.8 


1.161e-08 


63.5 


2.589e-09 


63.5 


4.788e-09 


63.8 


7r/240 


6.769e-09 


64.2 


1.785e-10 


65.0 


4.238e-ll 


61.1 


7.291e-ll 


65.7 


7r/480 


1.067e-10 


63.4 


1.750e-12 


102.0 


4.397e-12 


9.6 


1.558e-12 


46.8 



Table 2: Quadratic convergence of the maximum norm of the vector a, for the Kepler 
problem, by using the EHBVM(12,3) method, when imposing only the angular momentum 
conservation and both angular momentum and LRL vector conservation (a^^''). 



h 


al^^ ratio 


ajf^ ratio 


7r/30 

7r/60 

7r/120 

77/240 

7r/480 


4.530e-3 
1.155e-3 3.92 
2.902e-4 3.98 
7.265e-5 3.99 
1.837e-5 3.95 


1.246e-2 
3.195e-3 3.90 
8.040e-4 3.97 
2.013e-4 3.99 
5.055e-5 3.98 



• concerning the conservation of the angular momentum flM|) . all methods conserve this 
invariant, except HBVM(12,3), which is only energy-conserving. However, the error 
appears to be bounded, as is shown in Figure [3l 

• concerning the conservation of the LRL vector (1551) . all methods exhibit a drift, except 
EHBVM(12,3) when this invariant is required to be conserved. In particular: the drift 
of the 3-stage Gauss method and HBVM(12,3) method is practically the same, and 
larger than that shown by the EHBVM(12,3) method for which only the invariants 
(133!) and are imposed to be conserved. This, in turn, agrees with the analysis in 
[15] . As expected, no drift is observed, when the conservation of all the invariants is 
required for the EHBVM(12,3) method. 



5 Conclusions 

In this paper, we have used the technique of discrete line integrals introduced by lavernaro 
[23] to define an extension of the energy-conserving methods named HBVMs, in order to cope 
with the conservation of multiple invariants for Hamiltonian problems. This has resulted in 
an "enhanced" version of HBVMs, which we have, therefore, called EHBVMs. The analysis 
of such methods has been carried out, proving that the original order of HBVMs is retained 
by the new methods. At last, a few numerical tests clearly confirm the theoretical findings. 
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Figure 1: Components of the vector a. for the fully conservative EHBVM(12,3) method, 
h = 7r/30. 
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